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It is well established that the response of a black hole to a generic perturbation is characterized 
by a spectrum of damped resonances, called quasinormal modes; and that, in the limit of large 
angular momentum (I S> 1), the quasinormal mode frequency spectrum is related to the properties 
of unstable null orbits. In this paper we develop an expansion method to explore the link. We obtain 
new closed-form approximations for the lightly-damped part of the spectrum in the large-Z regime. 
We confirm that, at leading order in I, the resonance frequency is linked to the orbital frequency, 
and the resonance damping to the Lyapunov exponent, of the relevant null orbit. We go somewhat 
further than previous studies to establish (i) a spin-dependent correction to the frequency at order 
,— I l/l for equatorial (m = ±Z) modes, and (ii) a new result for polar modes (m = 0). We validate the 

approach by testing the closed-form approximations against frequencies obtained numerically with 
Leaver's method. 

Q-i I. INTRODUCTION 

U 

Quasinormal modes (QNMs) are damped resonances which play a key role in black hole dynamics [1-4]. For 
example, after the merger of a pair of black holes, the composite system 'rings down' through gravitational 
wave emission until it settles down into a quiescent axisymmetric Kerr phase. The gravitational wave signal 
f-~^ from ringdown is dominated by the least-damped QNM resonances. Chandrasckhar [5] was moved to compare 

the ringdown signal to the dying pure notes sounded by a bell. The analogy is apt; when a bell is struck, 
the dying tones depend only on the properties of the bell, rather than the hammer used to strike it, which 
affects only the relative degree of excitation of overtones. Likewise, the QNM spectrum depends only on the 
underlying properties of the black hole, rather than the complicated details of any initial perturbation. 

Physically, a black hole quasinormal mode is a decaying resonance which satisfies a pair of causally- motivated 
boundary conditions (typically being ingoing at the event horizon, and outgoing at spatial infinity). Mathe- 
matically, 'quasinormal ringing' emerges from a sum of residues of poles of the Green function in the complex 
frequency domain [6] . Each pole corresponds to a quasinormal mode of a single complex frequency. The real 
part of the frequency corresponds to the oscillation rate and the (negative) imaginary part corresponds to the 
damping rate. 

In the Kerr spacetime QNMs are labelled by three indices: multipole I, azimuthal number m and overtone 
number n > 0. The QNM spectrum depends only on the properties of the field (i.e. spin s) and the underlying 
black hole geometry (i.e. mass M, charge Q, angular momentum J = aM). In this paper, we concern ourselves 
with the most astrophysically-relevant case: the uncharged non-extremal Kerr black hole (0 < J < M 2 and 
Q = 0). 

A wide range of methods have been developed for determining QNMs. A non-exhaustive list [7] includes 
(i) time domain methods [8, 9]; (ii) direct integration in the frequency domain [10]; (iii) inverse potential 
methods [11, 12]; (iv) WKB methods [13-19]; (v) phase-integral methods [20-22]; (vi) continued fraction 
method [23, 24]; (vii) a semi-classical expansion method [25]; and (viii) an asymptotic iteration method 
[26, 27]. Review articles [1-4] explore some of these methods in greater depth. 
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The quasinormal frequencies of the Kerr black hole were first computed by Detweiler [28], and later by 
Leaver [23], Seidel and Iyer [17], Kokkotas [18] and Onozawa [29] (among others). In essence, the problem 
of accurately determining the QNMs of the Kerr black hole numerically has been 'solved' for a quarter of a 
century; the method introduced by Leaver [23] is fast, reliable and accurate, and with certain modifications 
[24], robust at large overtones n and angular momenta l,m. However, a numerical method provides little 
by way of physical insight. In contrast, approximation methods aim to provide some insight, perhaps at 
the expense of accuracy and applicability. For example, a body of work on the WKB method [13-17, 19] 
in spherically-symmetric spacetimes 'explains' the relationship between the low-overtone frequencies and the 
shape of the peak of the radial potential [15]. The low-overtone QNM spectrum of the Kerr black hole can 
be obtained via the WKB method [17, 18], although the expansion formulae are somewhat complicated. 

It has been known for many years that quasinormal modes are intimately linked to the existence and 
properties of unstable photon orbits [12, 30]. A quarter-century ago, Mashhoon [31] examined an aggregate 
of test null rays in orbit in the equatorial plane of a Kerr-Newman (Q > 0) black hole near the co-rotating 
circular orbit. Mashhoon's analysis suggested that, in the eikonal regime (I 3> 1), the maximally-corotating 
mode (m = I) has the frequency 

w{r'«*« + -*A + (n + l/2) (1) 

where w+ is the Kepler frequency for null rays on the unstable orbit and A+ is the decay rate of the unstable 
orbit (subsequently identified as the Lyapunov exponent [32-34]). These orbital quantities were plotted as a 
function of a in [35]; see also our Fig. 2. The 'geometric-optics' approach of Mashhoon is both appealing and 
insightful. For example, it neatly 'explains' why the maximally corotating modes (to — I) are faster-oscillation 
and slower-decaying than the other modes (to < I). However, Mashoon's original analysis was only 'indirect', 
in the sense that the Teukolsky equation, which describes massless perturbations of the Kerr spacetime, was 
not solved directly. So, for instance, it is not clear from (1) what effect the spin of the perturbing field has 
on the frequency spectrum. 

A recent work by Hod [36] demonstrated that a bridge may be built between "geometric-optics" intuition 
and "Teukolsky equation" exactitude. By analysizing the Teukolsky equation in the extremal limit (a — » M), 
Hod showed that (1) is indeed good approximation. 

The aim of this paper to show that improved approximations to QNM frequencies may be obtained by 
directly analysing the Teukolsky equation, without recourse to WKB methods. We go further than other 
studies (such as Mashhoon [31] and Hod [36]) to obtain frequency expansions in closed form which include the 
leading-order effect of the spin of the perturbing field, and which are valid outside the extremal regime (i.e. for 
all a). As well as considering the equatorial modes (to = ±1) we also analyse the 'polar' (m = 0) modes (see 
Fig. 1 for illustration) . We confirm that the expansions provide the correct fits to the numerically-determined 
spectrum, in the limit / 3> 1. 

To obtain our new analytic results, we apply a recently- introduced expansion method [25]. Previously, this 
method has been applied to a range of spherically-symmetric spacetimes [25], and to a simple cylindrically- 
symmetric system (the 'draining vortex') [37]. In this paper, we extend the method to treat an axisymmetric 
system for the first time. The key idea behind our approach is that the properties of modes of high angular 
momentum I >• 1 are related to the properties of the unstable photon orbits of the spacetime (see Fig. 1 
and [38]). This idea is not new [30]; the novel ingredient is a 'semi-classical' ansatz enables us build a 
bridge between the classical 'geodesic' picture of a battery of null geodesies near the photon orbit [31], and 
a perturbative analysis at the level of the Teukolsky equation. This ansatz allows the QNM frequencies and 
wavefunction to be expanded in inverse powers of I. 

The remainder of this paper is organized as follows: in Sec. II we describe the Kerr spacetime, null geodesies, 
and the Teukolsky equation; in Sec. Ill we develop the expansion method and obtain the key analytic results; 
in Sec. IV we validate the approach by testing the expansion against numerical results. We conclude in Sec. V 




FIG. 1: Critical null geodesies for a — 0.8A1. These plots show the special null geodesies which impinge from infinity 
and asymptote to the unstable photon orbits. The 2D plot (left) shows co- and counter-rotating orbits in the equatorial 
plane (with the spin axis pointing out of the page). The 3D plot (right) shows a polar orbit, with the spin axis of the 
BH pointing up the page. 



with a brief discussion. 



II. SPACETIME, GEODESICS AND WAVES 

In this section we lay some groundwork, and give geometrical motivation for the QNM analysis which 
follows in later sections. In Sec. II A we review some properties of the Kerr spacetime; in Sec. II B we consider 
the null geodesies of the spacetime, and the unstable photon orbits [38]; and in Sec. II C we recap the theory 
of weak-field perturbations of Kerr. 

A. The Kerr Spacetime 

The Kerr metric [39] describes the spacetime around a rotating uncharged black hole in vacuum (for details 
see e.g. [5, 40, 41]). In the Boyer-Lindquist coordinate system {t,r, 9, <j>} the line element takes the form 



2Ma 2 r sin 2 8 /p 2 )d(t) 2 , 

(2) 



ds 2 = -(1 - 2Mr/p 2 )dt 2 - (4aMr sin 2 0/p 2 )dtd(f> + (p 2 /A)dr 2 + p 2 d8 2 + sin 2 8(r 2 + a 2 

where A = r 2 — 2Mr + a 2 and p 2 = r 2 + a 2 cos 2 8, and where M and aM are the mass and angular momentum 
of the black hole. The spacetime has an ergosphere within r < 2M and horizons at r — r^± = M ±y/M 2 — a 2 . 

B. Null Geodesies and Orbits 



Killing vector fields (satisfying £( a; &) = 0) associated with time translation £?U and rotation around the 
symmetry axis £?,-, give rise to two constants of motion along a geodesic with tangent vector u a : an 'energy' 
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FIG. 2: Parameters for Unstable Photon Orbits around a Kerr black hole as a function of a — J/M . The left plot 
shows the orbital frequency with respect to coordinate time t (time measured at spatial infinity). The right plot shows 
the Lyapunov exponent, which controls the decay rate of QNMs. The lines are for the special cases of prograde and 
retrograde equatorial orbits (see also Fig. 1 & 2 in [35]) and the polar orbit, discussed in Sec. II B. General orbits 
would correspond to lines between these limits. 

E = —^f t \U a and an 'azimuthal angular momentum' L z = tf^Ua- The existence of a Killing tensor K a b 
satisfying Ku^-a = implies a further constant of motion, the 'Carter constant' [42] defined as Q — K a t,u a u b + 
L 2 — E 2 . The Carter constant governs the motion of geodesies in the latitudinal direction. The equations of 
motion for null rays can be written in first order form [5, 38, 40, 43] as 

(Y?E - 2MarL z ) /A, 
IMarE + {p 2 - 2Mr) 



P H 
P 2 j> 



(p 2 r)< 



P 



E 2 r 4 + {a l E 
L 2 



sin 2 

„2 



/A 
2M [{aE - L z 



Q]r-a 2 Q, 



= Q- 



sm 2 9 



-a 2 E 2 



(3) 
(4) 
(5) 
(6) 



where S 2 = (r 2 + a 2 ) 2 — a 2 A sin 2 9 and an overdot denotes differentiation with respect to an affine parameter. 
Note that the r and 9 equations may be decoupled into ordinary differential equations by using the 'Mino 
time' [43, 44] parameter v defined by dv = d\/p 2 . For our purposes this is not necessary. 

In the following section we consider three sets of constant-r null orbits: prograde and retrograde equatorial 
orbits with Q = and polar orbits with L z = 0. These orbits are illustrated in Fig. 1. The cases are henceforth 
distinguished using subscripts +, — and o, respectively; for example, we will denote the orbital radii by f + , 
f_ and f and the angular frequencies (with respect to coordinate time t) by u + , w_ and uj . Figure 2 shows 
the orbital frequency and Lyapunov exponent for these special orbits as a function of black hole rotation a. 

1. Equatorial orbits 



Consider first null geodesies in the equatorial plane, with Q — 0. The ratio of azimuthal angular momentum 
to energy defines an impact parameter (see Fig. 1), 



b = L z /E. 
Null geodesies in the equatorial plane are governed by the orbital equation 



b 2 



i 



2M(b - a) 



(7) 



(8) 



Perpetual photon orbits are possible if conditions r — and f — can be satisfied simultaneously. It is 
straightforward to show that both prograde (+) and retrograde (— ) orbits exist, at orbital radii r±, 

r± = 2M [1 + cos (| cos _1 (Ta/M))] , (9) 

with critical impact parameters b± given by 

b ± = ±3 y/Mf ± - a. (10) 

Note that we define the critical impact parameter b so that it is negative for the retrograde trajectory. The 
orbital frequency lj± — cf>/t of the circular orbit is obtained from the ratio of the <j> and t equations (3 and 4) , 
and it is straightforward to show that it is simply 

1 = MV 2 

UJ± \b±\ (f ± ) 3 / 2 ±aM 1 /2- { > 



In the critical case, the orbital equation (8) can be factorized into 

b± \ r J \ ' r J 



o-gu-±Y(^i m 



In [31, 32, 36] it was established that the decay rate for the equatorial circular orbits is controlled by 
Lyapunov exponent A (see Eq. (1)). For the equatorial orbits [32], it has been shown that the Lypanov 
exponent is 

A ± = l -' X± (13) 

\b ± \^l-x\ 

where x± is the ratio 

x± = a/b± (14) 

Note that 6_ < and hence £_ < 0. 

In the Schwarzschild limit (a — > 0) we recover r± = 3M, b± = \/27M and u>± = \± = l/(\/27M). In the 
extremal limit (a ->• M) we have r_ = AM, fe_ = -1M => cj_ = 1/(7M), A_ = 3/(7\/5M) and f+ = M, 
6 + = 2M =^> w + = 1/(2M) and A + = 0. Note that the decay rate of the co-rotating mode tends to zero in 
the extremal limit (as many authors have previously observed). 

2. Polar orbits 

Consider next the special case of trajectories with zero azimuthal angular momentum (L z = 0), which we 
will refer to as 'polar' geodesies. For polar geodesies we may define the natural impact parameter 

b 2 = Q/E 2 + a 2 (15) 

which corresponds to the perpendicular distance between the geodesic and the rotation axis (measured at 
spatial infinity). Note that polar geodesies are parallel to the rotation axis at spatial infinity, as shown in 
Fig. 1. 

Null polar geodesies are governed by the equation 

(£-y f) 2 ee R(r) =r 4 + (2a 2 - b 2 )r 2 + 2Mb 2 r + a 2 (a 2 - b 2 ) (16) 



An unstable orbit arises at radius r = f where R(f ) = and -R'(f ) — 0, given by 



f = M + 2v/M 2 - a 2 /3 



1 _, / M(M 2 -a 2 ) 
- cos ' 



3 V(^ 2 -a 2 /3) 3/2 

(see e.g. Eq. (14) in [38]). The critical impact parameter is given by 



(17) 



« = (3f °-.f K 1 + ° 2) . (18) 



and in the critical case b = b Ql the radial equation (16) can be factorized into 

R(r) = (r - f ) 2 (r 2 + 2f c r - a 2 (b 2 - a 2 )/f 2 ) . (19) 

The period T for a null geodesic on the unstable polar orbit at r = f can be found from integrating the 
ratio of the t and z — cos ((9) equations (see Appendix A for details), 

r +1 dt r +1 n+ a 2 z 2 /d 2 \ 1/2 

T = 2 / — dz = 2d — f± dz = Ad ellipE^o/do). (20) 



\ dz 7-i \ 1 — z 

where we have defined d = ^Jb 2 a — a 2 = y/Q/E, and ellipE(-) is the complete elliptic integral of the second 
kind [47] (a note of caution: ellipE(fc) is implemented in Maple as EllipticE(fc), but in Mathematica as 
EllipticE(fc 2 )). Hence the orbital frequency is 

7T 

W ° = 2d ellipE(m/d )' ^' 

We expect the imaginary part of the QNM of the polar modes (L z = 0) to be related to the Lyapunov 
exponent A for the polar null orbits. In Appendix A it is shown to be 



fo clhpK^o/y/l+^j / a 2 rf 2xV 2 

A ° (do) 2 yiT^ 2 cllipE(z a;o ) V ft J ' [Z2) 

where ellipK(-) is the complete elliptic integral of the first kind and we have defined the ratio 

x = a J do- (23) 

In the extremal limit, a -> M, we have f /M = 1 + y/2 w 2.4142 and b Q /M w 4.8284 and Mlu ~ 0.20937 
and MX w 0.162006. Unlike the corotating modes, the m = modes are still significantly damped in this 
limit (as is shown in Fig. 2). 

3. General orbits 

Above we considered two limiting cases: equatorial orbits restricted to 9 — ir/2, and polar orbits that 
explore the full range of polar angles, < 9 < it. A general orbit will explore a more limited range of polar 
angles, tt/2 — A9 < 9 < tt/2 + A9. The angular range A9 is a function of the ratio L z /y/Q, with limiting 
cases A9(L z /y/Q — $• oo) = and A9{L z j\j ' Q — \ 0) = tt/2. It was shown in e.g. [38] that perpetual orbits 
exist for all ratios L z / ' \[Q. A general analysis is left for future work. 

C. Perturbations of the Kerr Spacetime 

To analyse gravitational dynamics around a rotating black hole, one may examine perturbations of the 
metric tensor via Einstein's equations linearized around the Kerr background (assuming the perturbations are 
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"small"). This approach yields a set of coupled partial differential equations, which are difficult to analyze. 
An alternative approach is to consider perturbations of the Weyl and Ricci scalars using the Newman- Penrose 
formalism. This was the approach taken by Tcukolsky [45, 46], who decoupled and separated the equations 
to reduce them to a single 'master equation' for the radial part of the perturbation. 

A 'master equation' describes masslcss perturbations of scalar (|s| = 0), spinor (|s| = 1/2), electromagnetic 
(\s\ = 1) and gravitational (\s\ = 2) types. In vacuum, the radial equation is 

d ( +1 dR\ i ( K 2 - 2is(r - M)K 



A ~'* r + VJ + [ A- '-+**»-*)* = <> (24) 

(Eq. (4.8) in [46]) where K = (r 2 + a 2 )ui — am. The angular separation constant A = A[ m — 2maoj + a 2 ui 2 is 
determined from imposing physical boundary conditions at the poles (6 = 0, 7r) on the angular equation 



d 

dz 



*-< 



a u> z — 2amojsz + s + Ai m — 



(rn + sz) 
\-z 2 



5 = 0, where z — cos 8, (25) 



whose solutions are oblate spin-weighted spheroidal harmonics S = s Si m (au>; z) [47]. An expansion of the 
angular eigenvalue in powers of auj is given in [48, 49]. 

Let us make the substitution R(r) = r _1 A _s u(r) to write the radial equation (24) as 



d 2 u 2s(r - M) du 
dr 2 r 2 dr* 



n 2 - 2ls(r ~ 



is{r-M)n /A /' 2Ms Aicus\ 



u = (26) 



where 



and 



f(r) = A/r 2 = 1 - 2M/r + a 2 /r 2 (27) 



0(r) = K(r)/r 2 = (1 + a 2 /r 2 )uj - am/r 2 . (28) 

Here we have defined a tortoise coordinate r* via 

^ = r 1 . (29) 

dr 

Note that this definition differs from, e.g., Tcukolsky [46], who defines instead ^r = (r 2 + a 2 )/ A. 

Quasinormal modes are the complex-frequency modes that are purely ingoing at the (outer) horizon (r = 
rh+), and purely outgoing at spatial infinity, satisfying boundary conditions 

u (r) ~ { exp ( _i [ w _ am l r l+V*)' r * "»■ -°°. ( 30) 

I exp(+iwr sf ), r* — > +oo. 

III. QNM EXPANSION METHOD 

In this section we extend and develop the QNM expansion method first described in [25]. The method has 
three steps: (i) introduce an ansatz for the wavefunction, inspired by the null geodesies near the perpetual 
orbit, (ii) expand the wavefunction and frequency in inverse powers of in (or I), (iii) impose a continuity 
condition at the perpetual orbit radius (i.e. r = f± or r = f ) to recover the expansion coefficients. 



A. Equatorial modes of scalar field 



Let us illustrate the method by considering first the maximally corotating modes, m = I, of the scalar wave 
equation (,s = 0). We seek solutions to (26) that satisfy the QNM boundary conditions (30). Step (i) is to 
propose an ansatz of the form 



(31) 



u(r) = exp i / {3(r)dr* I v (r) 



We demand that /3(r) has three properties. First, f3(r) should have a root at the circular orbit r = f +1 so 
that /3 changes sign here, and the solution passes smoothly from outgoing behaviour at infinity to ingoing 
behaviour at the horizon. Second, 1 should be such that, upon insertion into Eq. (26), the resulting equation 
has an overall factor of /. That is, 



-(3 2 + O 2 ex /, 



(32) 



Third, /? should ensure that the QNM boundary conditions are satisfied, provided v is regular at both horizon 
and infinity. 

The third property suggests that ft has an overall factor of f2(r) (defined in Eq. (28)), and the first property 
suggests that we take inspiration from the factorized form of the orbital equation (12) to write 



fi(r) = O(r) 1 



a(b — a) 



1 



1 



2n 



1/2 



(33) 



where b = 6+, defined in (10). It is straightforward to confirm that the second property (32) is satisfied, 

2 



-p 2 + n 2 = f(r)n 2i ^^- 



l- 



a(b — a) 



(34) 



and that the QNM boundary conditions are satisfied if v(r) is regular at the horizon and infinity. Upon 
substitution of the ansatz (31) into (26), and after dividing through by /, we reach the new radial equation 



d / dv 

dr \ dr 



2i(3 



dv 
dr 



n 



2 (b-a) 2 



1- 



a(b — a) 



■*"?-7 



v = 0. 



(35) 



Step (ii) is to expand the wavefunction, frequency and angular eigenvalue in inverse powers of m = I, in 
the following manner, 



buj = txj-im + wo + vo\m +... 

v(r) = exp (So(r) + m~ 1 Si(r) + . . .) 

A = A_ 2 ra 2 + A_ito + A + . . . 



(36) 

(37) 
(38) 



Here, {ti7_ 1 , zuq, . . .} and {A_ 2 > A_i, . . .} are coefficients to be determined, and {So(r), S\(r), . . .} are regular 
functions of r. Inserting (36, 37, 38) into (26) and grouping together like powers of m leads to a system of 
equations, 



0(m 2 ) 



(b-a) 2 
b 2 



1 



G7_l 



a b 



1 - 



a(b — a) 



-A_ 2 =0 



0(mi): 2^f^ + ^-^i + 2(l + ^) 
O(m ) : ... 



2X _ (b-af 

WO TTT^ = 



b 2 r 2 



(39) 

(40) 
(41) 



Here /3 (r) = (zu-i/b)(l - a(b - a)/r 2 ) _1 (l - f+/r)(l + 2f + /r) 1/2 and 6 = b + given in Eq. (10). A major 
difference with the Schwarzschild case is that the angular expansion coefficients A/, is also dependent on 
frequency, through acu. This challenge is not insurmountable. First we note that A_2 = (1 — anj_i/6) 2 . It is 
trivial to then show that the choice w—\ — 1 satisfies Eq. (39). To find higher-order coefficients, we may use 
the series expansion of the eigenvalue given in [48, 49]. We find 

A_ 2 = (l-x) 2 (42) 

A_i = 5i - 2x(l - x)m (43) 

A = -2s 2 xS 2 + —S 3 + x—±m (44) 

4 ax 

where x = a/b and 

Sx = 1 - §x 2 - §x 4 - 1/ + 0(x s ) = (1 - x 2 ) 1/2 (45) 

S 2 = l + .T + x 2 +x 3 + x 4 + 2; 5 +e>(.T 6 ) = (1 - x)" 1 (46) 

S 3 = l + x 2 +x 4 + e>(x 6 ) = (1-x 2 )- 1 (47) 

Though Seidel [48] provides an expansion to sixth order in au, this provides only the first few terms of a 
(presumably) infinite series. Above, we used intuition to 'guess' a closed form from the first few terms of the 
series expansion (the results to the right of the = symbol) . Mathematica can be used to obtain further terms 
in the expansion of the spheroidal eigenvalue in the scalar case (s = 0); the higher terms are consistent with 
the assumption. In Sec. IV we test the resulting frequencies against numerical results and find the expected 
agreement. 

Next we demand that the radial function So(r), which features in the expansion of the wavefunction (37), 
is continuous at r = r+ (where /3(f+) = 0). After inserting result (43) into (40) and evaluating at r = r+ and 
rearranging, we obtain the coefficient 

-.-^,(1-0 (48) 

The function So may then be determined by substituting (48) back into (40) and rearranging to find S , then 
integrating. 

In principle, we may continue in this fashion to determine the higher-order coefficients vd\ , w-i , . . . and phase 
functions Si(r), ^(r), . . ., with the help of a symbolic algebra package. In practice it is a challenge to continue 
the expansion beyond 0(m~ 2 ), due to the frequency dependence of the angular eigenvalue. Results to this 
order are given in Eq. (51). 

To look for higher modes (n > 0) one must first modify the ansatz for the wavefunction (37), as shown in 
[25]. This technique follows through to the axisymmetric case without additional difficulties, although the 
calculation is not pursued here. 



B. Equatorial modes of higher spin 

Let us now show how one may generalise the previous analysis to treat fields of higher spin (s^0). As in 
Sec. Ill A, we start with Eq. (26) and once again propose an ansatz of the form (31). As in Sec. Ill A, we 
demand that /3 obeys three conditions, but now the second condition (32) becomes 

/ is(r-M)\ 2 / is{r-M)\ 2 

-\P+ 2 + n- y 2 ' ) ex/, (49) 
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The conclusion is that 

t> + l3{r ; 2 M) = (n - " (r - M) ) (i - °^) " (i - ^) (i + ^) V2 . (50) 

It is straightforward to verify that /3 satisfies the QNM boundary conditions at the horizon and at infinity. 

The method proceeds through step (ii) and (iii) as in the scalar case. The method was automated using a 
symbolic algebra package. We now present the key results for the equatorial modes. 

C. Equatorial Modes: Key Results 

Applying the method described above, we find that the fundamental (n = 0) maximally co-rotating (m = I) 
frequency w„™^~ has the following 'eikonal' expansion : 

b + 43 = m + ^=^(1 - + 21fi ^ ^ )2 [(7 + 44, + 127x 2 ) - 72 s 2 (l + xf] + 0(m- 2 ) (51) 

where here x = x + = a/b + and the critical impact parameter b + was defined in (10). We established in 
Eq. (11) that 1/6+ is equal to the Kepler orbital frequency w + for null rays in corotating circular orbit. 
Furthermore, the leading-order imaginary component is equal to one-half of the Lyapunov exponent given in 
(13) (and Refs. [31, 32]). 

To leading order, our result is consistent with the outcome of the geodesic analysis of Mashhoon [31], 
Cardoso et al. [32] and Hod ([36], Eq. 3 and 4), and with the wave-equation analysis of Hod in the extremal 
regime. Note however that result (51) goes two steps further. Firstly, it includes a spin-independent correction 
to the real part of frequency at order O(m ), which has not been previously obtained. Secondly, it also includes 
the spin-dependent correction at order 0(m~ l ) for the first time. In Sec. IV we compare result (51) against 
numerically-determined frequencies, to verify that both new terms are correct. 

The expansion of the fundamental (n = 0) maximally counter-rotating (m = —I) frequency lo^Zq is 

b- 4 m n =o l) = ™- ^=S(! - + 216 ( j l ( 1 ^y [(^ + 44^ + 127x 2 ) - 72,s 2 (l + X f] + 0(m^) (52) 

where here x = X- = a/6_. Note that the counter-rotating impact parameter b- given by (10) is negative, as 
is x (for positive a). Result (52) may be obtained from result (51) via the simultaneous replacements b + — > 6_, 

-i (m—l) (m=— I) 

x + — > x^ , m — > —m and w„ =0 — ^ — w n=0 

The QNM spectrum has the following symmetry property: 

w lmn = ^^l-ra.n- (53) 

Therefore frequencies with negative real part may be obtained by applying (53) to Eqs. (52) and (51). 

D. QNMs of polar modes (m = 0) 

In Sec. II B 2 we considered 'polar' geodesies, which have zero azimuthal angular momentum L z = and 
pass through the north and south poles 9 — 0, tt. Here we apply the expansion method to the Teukolsky 
equation to find the corresponding m — QNM frequencies. 

For brevity, let us consider scalar waves (s = 0) here, as it turns out that the leading order terms (at I 1 
and 1°) are not spin-dependent. The appropriate ansatz for the polar modes is again of the form (31), with 
condition (ii) (Eq. (32)) suggesting the choice 

f = rt-f { r)^ = ^, (54) 
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where b Q is the critical impact parameter (18) and R(r) is the quartic defined in Eq. (16), with factorization 
given by Eq. (19). 

Inserting the ansatz (54) into (26) leads to 



C/V )' + 2»0(r)v / + 



•.«' + ^4 



v = 0. 



Next we expand in inverse powers of L = I + 1/2, 

b uj = m_\L + zdq + W\h + . . . 
v(r) = exp(5 (r) + 5i(r)L- 1 + ...) 

A = A- 2 L 2 + A-iL + A + . . . 

Now group terms in (55) order-by-order in L, 



n7 2 _ 1 - A_ 



2in7_i (r 2 + 2rf - a 2 d 2 /f 2 ) 



(r-f o )S' (r) + 



IW—XVJq 



■i$> 



A_i 



= 

= 
= 



(55) 



(56) 

(57) 
(58) 



(59) 

(60) 
(61) 



Once again, there remains an obstacle to progress: the frequency-dependence of the angular eigenvalue A. 
The angular eigenvalue of the scalar field (s = 0) for the polar mode A; m=0 has the following expansion: 



A 



l,m=0 



L 2 



c 2 _ l 

2 4 



4c 2 5 c 4 - 8c 2 5c 8 - 160c 6 + 2256c 4 - 1024c 2 



32L 2 



64 L 4 



8192L 6 



■*°ljj>)> (62) 



where c~auj and L = I + 1/2. This result is given in Ref. [50], Theorem 10. For higher spin fields s/Owc 
may use the expansions given in [48, 49]: 



A = L 2 



1 c 2 1 c 4 / c» 

1+ 2tf + 322? +0 + °U 



, ^ 1/A (I 2 \c 2 /5 3s 2 \ c /-5 5s 2 \ c 6 /c 8 

This expansion allows us to write A_2 and A_i as power series in aro_i, where 

A-2 = l + ^aw_ 1 ) 2 + ^(a W _ 1 ) 4 + + ^(au7_ 1 ) 8 + O((au7_ 1 ) 10 ), 



A_i = (au7_i)(aw ) 



and 



1 + \(aw^) 2 + + ^(aro.!) 6 + 0((ao7_i) 8 ) 



a = a/6 Q - 



(63) 

(64) 
(65) 

(66) 



Hence to find W-i we must solve a non- linear equation. This is straightforward to do iteratively. We may 
obtain a converging sequence of estimates {w_i = 1, tuLx, ro_i, • ■ •} using 



To find wo we impose a continuity condition on So(r) a ^ r — ^o, and solve to obtain 
-i(3r 2 - o 2 ^/f2)V2 



G7 



2/Jo 



1 - — I 1 + -(cra7_i) 2 + H (cra7_i) 6 + . . 

2 V 8 V ' 1024 V l> 



1/2 



-i -1 



(67) 



(68) 
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Note that vj-\ is purely real and tuq is purely imaginary. 

It is unfortunate that closed-form results arc not so easily obtained in the polar (to = 0) case. However, we 
may compare Eq. (67) with Eq. (21), and Eq. (68) with Eq. (22), by evaluating numerically for a given rotation 
parameter a. It turns out that the numerical values are in precise agreement. For example, for a = 0.8M 
we have f Q w 2.67062M, b w 4.98488M and d w 4.92027M. Iteration of Eq. (67) leads to ra7_i/& sa 
0.20191308580895 795M~\ which should be compared with uj « 0.20191308580895 903M- 1 from Eq. (21), 
i.e. agreement to 14 significant figures. Equation (68) leads to wo/b = — 0.0893785765669 3478zM _1 , which 
should be compared with -i\ Q /2 = -0.0893785765669 4054iM" 1 from Eq. (22), i.e. agreement to 12 signifi- 
cant figures. 

We conclude that the polar (to = 0) modes have the following frequency expansion: 

"£=? = U {l + l/2) - l\ Q /2 + 0{l7l- 1 ) 

= (w_i/6 e )(i + 1/2) + (tn /6 o )(2n + l)+O(m- 1 ) (69) 

where {uj , A } are given in Eq. (21) and (22), and {w-i,zuo} are given in Eq. (67) and (68). In Sec. IV we 
compare this estimate against numerically-determined frequencies. 



IV. VALIDATION 

In this section we test the key results (51), (52) and (69) by comparing with numerically-determined 
frequencies. A fast and accurate numerical method for determining QNM frequencies for the Kerr black hole 
was introduced by Leaver [23] many years ago. We have implemented our own version of the continued-fraction 
algorithm that Leaver described. Other implementations are available in the public domain [4, 51]. 

1. Equatorial Modes 

Table I shows the fundamental (n = 0) frequencies for scalar-field (s = 0) modes with \m\ — I — 2, 4, 6, 8, 10 
at a = 0.8M. For each I, the upper row gives the numerically-determined frequency The lower rows give 
the approximations at 0th (to 1 ), 1st (to ) and 2nd (m 1 ) orders, obtained from Eq. (51) and (52). Previous 
approximations [31, 32, 36] only supplied the real part to order m 1 (0th), and the imaginary part to order 
m° (1st). The table shows that including the higher-order corrections to the real part improves the accuracy 
of the estimate substantially. In fact, even at relatively low I, Eq. (51) and (52) are surprisingly accurate 
estimates of the scalar-field frequencies. For example, at I = 2 the estimate of the real part is accurate to 
0.3% (co-rotating) and 0.06% (counter-rotating). The estimate of the imaginary part is accurate to 1.1% 
(co- rotating) and 0.6% (counter-rotating). 

Table II shows the fundamental (n = 0) frequencies for the gravitational field (|s| = 2) at a = 0.8M. Again, 
it is clear that including the higher-order corrections improves the accuracy of the estimate. However, the 
magnitude of the error in the estimate is significantly greater than for the scalar field, implying that field 
spin has a non- negligible effect outside the cikonal regime (i.e. at small or moderate /). For example, the real 
part of the I — 2 estimate is in error by 10.4% (co-rotating) and 2.2% (counter- rotating) , and by 6.6% and 
5.9% for the imaginary part. In the \ow-l regime the expansion is substantially less accurate than the WKB 
method [17, 18]. 

Figure 3 shows the 'error' (defined as difference between the estimate given in Sec. Ill and the numerically- 
determined frequency) as a function of to, on a log-log scale. It provides strong evidence that the estimates 
given in (51) and (52) are indeed correct to the stated order in m. The upper plots show the error in the 
real part of frequency for the scalar (left) and gravitational (right) cases, for co-rotating orbits. The middle 
plots show the same for the counter-rotating orbits. The data set marked "0th" shows the error using only 
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Prograde 


(m = l) 


Retrograde 


(m = -I) 


\m\ 




Re(Muj) 


Im(Mw) 


Re(Mw) 


Im(Mw) 


M = 2, 


exact 


0.70682338 


-0.08152026 


0.39573382 


-0.09428527 




2nd 


0.70892744 




0.39550921 






1st 


0.69841410 


-0.08061490 


0.39393550 


-0.09374776 




0th 


0.61779920 




0.30018774 




\m\ = 4, 


exact 


1.32089554 


-0.08088308 


0.69498056 


-0.09393068 




2nd 


1.32146997 




0.69491010 






1st 


1.31621330 


-0.08061490 


0.69412324 


-0.09374776 




0th 


1.23559839 




0.60037549 




\m\ = 6, 


exact 


1.93725163 


-0.08074051 


0.99486912 


-0.09383851 




2nd 


1.93751694 




0.99483556 






1st 


1.93401249 


-0.08061490 


0.99431099 


-0.09374776 




0th 


1.85339759 




0.90056323 




\m\ = 8, 


exact 


2.55428739 


-0.08068741 


1.29491167 


-0.09380178 




2nd 


2.55444002 




1.29489216 






1st 


2.55181169 


-0.08061490 


1.29449873 


-0.09374776 




0th 


2.47119678 




1.20075097 




|m| = 10, 


exact 


3.17161440 


-0.08066204 


1.59501395 


-0.09378355 




2nd 


3.17171355 




1.59500121 






1st 


3.16961088 


-0.08061490 


1.59468647 


-0.09374776 




0th 


3.08899598 




1.50093872 





TABLE I: Equatorial Modes: Scalar Field. For each \m\, the top row is the numerically-determined QNM frequency 
of the fundamental mode, for m = I (left) and m = —I (right). The lower rows give the estimates from Eq. (51) and 
(52) at orders m 1 (0th), m° (1st) and mT 1 . 







Prograde 


(m = l) 


Retrograde 


(m = -I) 


\m\ 




Re(Muj) 


Im(Mw) 


Re(Mu;) 


Im(Mw) 


M = 2, 


exact 


0.58601697 


-0.07562955 


0.30331342 


-0.08851224 




2nd 


0.52518088 




0.29659659 






1st 


0.69841410 


-0.08061490 


0.39393550 


-0.09374776 




0th 


0.61779920 




0.30018774 




\m\ = 4, 


exact 


1.24754701 


-0.07812549 


0.64854056 


-0.09258611 




2nd 


1.22959668 




0.64545379 






1st 


1.31621330 


-0.08061490 


0.69412324 


-0.09374776 




0th 


1.23559839 




0.60037549 




\m\ = 6, 


exact 


1.88474955 


-0.07927839 


0.96349484 


-0.09323669 




2nd 


1.87626808 




0.96186468 






1st 


1.93401249 


-0.08061490 


0.99431098 


-0.09374776 




0th 


1.85339759 




0.90056323 




M = 8, 


exact 


2.51341639 


-0.07979031 


1.27116161 


-0.09345939 




2nd 


2.50850338 




1.27016400 






1st 


2.55181169 


-0.08061490 


1.29449873 


-0.09374776 




0th 


2.47119678 




1.20075097 




\m\ = 10, 


exact 


3.13816080 


-0.08005731 


1.57589011 


-0.09356231 




2nd 


3.13496424 




1.57521869 






1st 


3.16961088 


-0.08061490 


1.59468647 


-0.09374776 




0th 


3.08899598 




1.50093872 





TABLE II: Equatorial Modes: Gravitational Field. As Table I but for the gravitational QNMs. 
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I 


Re(Moj) 


Im(Mw) 


1 = 2, 


scalar 


0.50712267 


-0.08966931 




gravitational 


0.40191735 


-0.08215627 




estimate 


0.50478271 


-0.08937858 


I = 4, 


scalar 


0.90990916 


-0.08947052 




gravitational 


0.85459105 


-0.08741415 




estimate 


0.90860889 


-0.08937858 


1 = 6, 


scalar 


1.31333515 


-0.08942290 




gravitational 


1.27557678 


-0.08848090 




estimate 


1.31243506 


-0.08937858 


1 = 8, 


scalar 


1.71694949 


-0.08940455 




gravitational 


1.68823375 


-0.08886365 




estimate 


1.71626123 


-0.08937858 


1 = 10, 


scalar 


2.12064455 


-0.08939561 




gravitational 


2.09746074 


-0.08904433 




estimate 


2.12008740 


-0.08937858 



TABLE III: Polar Modes. The table compares the numerically-determined values of the 'polar' (m = 0) modes (top 
two rows) with the lowest-order estimate (69). Note that at leading order the estimate is spin-independent. 

the order m 1 estimate. The data sets marked "1st" and "2nd" show the effect on error on including the 
order m° and order mT 1 corrections. The plots shows that the error scales as m° (0th), to -1 (1st) and 
m~ 2 (2nd) in the large-m regime (this may be inferred by examining the 'slope' of the respective data sets 
in the log-log plot, and confirming that it tends to 0, —1 and —2 in the large-^ limit, for 0th, 1st and 2nd 
approximations). Comparing the left and right plots, it is clear that the absolute error is significantly greater 
for the gravitational field than for the scalar field. 



2. Polar Modes 

Table III compares numerically-determined polar (m = 0) frequencies for scalar and gravitational fields, 
with the approximation of Eq. (69). The approximation is closer to the scalar QNM frequency than to the 
gravitational QNM frequency. In the lower plots of Fig. 3, the error in the real and imaginary part of frequency 
is shown. It is clear from the gradients of the two data sets that the next-order correction to (69) is at 0(l^ r ) 
for the real part, and 0(l~ 2 ) for the imaginary part. 

V. DISCUSSION AND CONCLUSIONS 



In this paper we have obtained new closed-form approximations for the fundamental (n = 0) Kerr QNM 
frequencies of the equatorial (m = ±Z) and polar (m = 0) modes of arbitrary spin, valid in the regime I ^$> \s\; 
see Eq. (51), (52) and (69). The result for the corotating modes improves upon previous approximations 
[31, 32, 36] in two regards, giving: (i) a spin-independent correction to the real part of the frequency at order 
m , (ii) a spin-dependent correction to the real part at order mT x . The result for the polar modes is (we 
believe) substantially new. We demonstrated that the approximations fit the numerical data in the asymptotic 
regime. We have sought to improve understanding of the link between the asymptotic structure of the QNM 
spectrum and the properties of the unstable photon orbits (see also [30-32, 36]). 

Let us conclude by suggesting some possible extensions to this work: (i) it is simple to extend the analysis 
to higher overtones (n > 0) following [25] ; (ii) it should be straightforward to extend the analysis to the Kerr- 
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Accuracy of Frequency Estimate for Co-rotating Equatorial Mode: Scalar Field 



Accuracy of Frequency Estimate for Co-rotating Equatorial Mode: Gravitational Field 



a = 0.8M ' 

lsl=0 * „ 

m - I (corotating) * 

Oth order estimate + 
1 st order estimate 
2nd order estimate * 

10 
azimuthal number m 

Accuracy of Frequency Estimate for Counter-rotating Equatorial Mode: Scalar Field 



a = 0.8M 
Isl = 

m - -I (counter- rotating) 

Oth order estimate + 
1 st order estimate 
2nd order estimate ■ 



azimuthal number Iml 
Accuracy of Frequency Estimate for Polar Mode: Scala 



0.0001 






a = 0.8M 




Isl = 




m - (polar 


1e-05 


- Bejo) * 
lm(ffij X 



a = 0.8M 

Isl -2 

m - I (corotating) 

Oth order estimate + 
1st order estimate 
2nd order estimate * 

10 
azimuthal number m 

Accuracy of Frequency Estimate for Counter- rotating Equatorial Mode: Gravitational Field 



a = 0.8M 

Isl -2 

m - -I (counter-rotating) 

Oth order estimate + 
1st order estimate 
2nd order estimate * 

10 
azimuthal number Iml 

Accuracy of Frequency Estimate for Polar Mode: Gravitational Field 



a = 


0.8M 


Isl 


= 2 


m 


= (polar 


Re(u 
lm(u 


) + 

) x 



angular momentum, 1+1/2 



angular momentum, 1+1/2 



FIG. 3: Accuracy of the QNM Frequency Estimates. These log-scale plots show the difference between the numerically- 
determined QNM frequency and the estimates obtained in this paper, as function of angular momentum I. The left 
plots show the scalar field, and the right plots show the gravitational field. The upper (m = I) and middle (m = — I) 
plots show the difference in the real part of the frequency, for approximations (51) and (52) truncated at orders I 
(Oth order), 1° (1st) and l^ 1 (2nd). At large I, the slopes of the data sets on the log-log scale approach 0, —1 and —2, 
as expected. The lower plots show the difference for the real and imaginary parts of the polar (m = 0) frequencies. 
The gradients tend to -1 [red, real] and -2 [blue, imaginary], implying that the approximation is accurate up to order 
0(L _1 ) for the real part and 0(L~ 2 ) for the imaginary part. 
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Newman (Q > 0) black hole (though here coupling between electromagnetic and gravitational perturbations 
may complicate the issue [52, 53]); (iii) here we have considered only the limiting cases of equatorial (|m| = I) 
and polar (m = 0) modes here; extension to general modes (m = /j,l, —1 < n < 1) remains an open challenge. 
We expect the real and imaginary parts to again be linked to the orbital frequency and Lyapunov exponent 
on the relevant null orbit (i.e. the general orbits considered in [38]), although we also expect some non-trivial 
correction at O(l ), such as appears in (51); (iv) The expansion method detailed here could also be applied to 
find the Regge poles of the Kerr black hole in the eikonal regime [54, 55]; (v) The expansion method can also 
be applied to compute QNM excitation factors [6, 56], which are needed to tackle the question 'how much 
QNM ringing is excited by a given perturbation?' We hope to pursue this calculation in the near future. 

Appendix A: Period and Lyapunov Exponent for Polar Orbit 



Here we derive expressions for the orbital frequency and Lyapunov exponent for the polar orbit. The period 
T Q with respect to the coordinate time is 



T = 2 



1 i 



dz, 



where, from Eq. (3) and (6), 



p 2 A 



i=p- 2 (^ +a 2 z 2 ) l/2 (1 _ z 2 ) l/ S 



and do = y/b — o? = \fQjE. The t expression can be simplified with the use of the identities 



r 2 - a 2 \ 1/2 
r = d ( ° 2 ) , M = f 



( ?l + a 2 
\ 3f o — a 2 



to obtain 



Hence the period is 



i = p- 2 (dl + a 2 z 2 ) 



f 1 (d 2 + a 2 z 2 ) 1/2 

r = 2/ l J* dZ 



-i (1-z 2 ) 

which evaluates to result (20). 

The Lyapunov exponent A Q can be found from taking an orbital average, 



Xo = iJj{l u '- 



of the expression [32] , 



A= t 



;-i 



ld 2 V r 



2 dr 2 



where V r = r 2 = p 4 (r - f a ) 2 {r 2 + 2rr — a 2 d 2 jf 2 ). Using V" = 2p 4 (3f 2 — a 2 d 2 /f 2 ) gives 

2(3f 2 - a 2 d 2 Jf 2 ) l l 2 f 1 dz 



A = 
which leads directly to result (22 



T 



i (d 2 +a 2 z 2 ) 1/2 (l -z 2 ) 1 ' 2 ' 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 
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